---
title: "R Notebook"
output: html_notebook
---

This is an [R Markdown](http://rmarkdown.rstudio.com) Notebook. When you execute code within the notebook, the results appear beneath the code. 

Try executing this chunk by clicking the *Run* button within the chunk or by placing your cursor inside it and pressing *Cmd+Shift+Enter*. 

#read in data and subset
```{r}
library(contrast)
test = read.table("female EAD active compounds_hm_all samples_norm2.csv", sep=",", header=TRUE)
nor = subset(test, test$Pop=="N")
can = subset(test, test$Pop=="C")

#remove canadian outliers
can = can[-c(56,93),]

```


#analyze norwegian data separately
```{r}

start <- 5
end<- 24

for (i in start:end) {
  
lm_fit <- glm(log(nor[,i]) ~ Sex*Genotype, data = nor)  

sum = summary(lm_fit)

males <- contrast(lm_fit, 
                     list(Sex = "M", Genotype = "AA"),
                     list(Sex = "M", Genotype = "BB"))

females <- contrast(lm_fit, 
                     list(Sex = "F", Genotype = "AA"),
                     list(Sex = "F", Genotype = "BB"))




results = cbind(names(nor[i]),"Norway",sum$coefficients[2,3],sum$coefficients[2,4],sum$coefficients[3,3],sum$coefficients[3,4], sum$coefficients[4,3],sum$coefficients[4,4], males$testStat,males$Pvalue, females$testStat, females$Pvalue )



  
if (i==start) {
  
        
  output<-as.data.frame(results)

  }
  
  
  if (i>start) {
    output= rbind(output,results)

  
  } 
}

colnames(output)<- c("Compound","Population", "Sex_t","Sex_p","Geno_t","Geno_p","GxS_t","GxS_p", "male_contrast_t", "male_contrast_p", "female_contrast_t", "female_contrast_p")


output$Sex_p = as.numeric(as.character(output$Sex_p))
output$GxS_p = as.numeric(as.character(output$GxS_p))
output$male_contrast_p = as.numeric(as.character(output$male_contrast_p))
output$female_contrast_p = as.numeric(as.character(output$female_contrast_p))



output$sex_fdr = p.adjust(output$Sex_p, method = "fdr", n = length(output$Sex_p))
output$GxS_fdr = p.adjust(output$GxS_p, method = "fdr", n = length(output$GxS_p))
output$male_contrast_fdr = p.adjust(output$male_contrast_p, method = "fdr", n = length(output$male_contrast_p))
output$female_contrast_fdr = p.adjust(output$female_contrast_p, method = "fdr", n = length(output$female_contrast_p))


output$Geno_p = as.numeric(as.character(output$Geno_p))
output$Geno_fdr = p.adjust(output$Geno_p, method = "fdr", n = length(output$Geno_p))
  
```




#analyze canadian data separately
```{r}

start <- 5
end<- 24

for (i in start:end) {
  
lm_fit <- glm(log(can[,i]) ~ Sex*Genotype, data = can)  

sum = summary(lm_fit)

males <- contrast(lm_fit, 
                     list(Sex = "M", Genotype = "AA"),
                     list(Sex = "M", Genotype = "BB"))

females <- contrast(lm_fit, 
                     list(Sex = "F", Genotype = "AA"),
                     list(Sex = "F", Genotype = "BB"))




results = cbind(names(can[i]),"Canada",sum$coefficients[2,3],sum$coefficients[2,4],sum$coefficients[3,3],sum$coefficients[3,4], sum$coefficients[4,3],sum$coefficients[4,4], males$testStat,males$Pvalue, females$testStat, females$Pvalue )



  
if (i==start) {
  
        
  output_c<-as.data.frame(results)

  }
  
  
  if (i>start) {
    output_c= rbind(output_c,results)

  
  } 
}




colnames(output_c)<- c("Compound","Population", "Sex_t","Sex_p","Geno_t","Geno_p","GxS_t","GxS_p", "male_contrast_t", "male_contrast_p", "female_contrast_t", "female_contrast_p")

output_c$Sex_p = as.numeric(as.character(output_c$Sex_p))
output_c$GxS_p = as.numeric(as.character(output_c$GxS_p))
output_c$male_contrast_p = as.numeric(as.character(output_c$male_contrast_p))
output_c$female_contrast_p = as.numeric(as.character(output_c$female_contrast_p))


output_c$sex_fdr = p.adjust(output_c$Sex_p, method = "fdr", n = length(output_c$Sex_p))
output_c$GxS_fdr = p.adjust(output_c$GxS_p, method = "fdr", n = length(output_c$GxS_p))
output_c$male_contrast_fdr = p.adjust(output_c$male_contrast_p, method = "fdr", n = length(output_c$male_contrast_p))
output_c$female_contrast_fdr = p.adjust(output_c$female_contrast_p, method = "fdr", n = length(output_c$female_contrast_p))


output_c$Geno_p = as.numeric(as.character(output_c$Geno_p))
output_c$Geno_fdr = p.adjust(output_c$Geno_p, method = "fdr", n = length(output_c$Geno_p))
  
```

#combine and output
```{r}
full = rbind(output, output_c) 

write.table(full, "model_output.csv", sep=",", row.names = FALSE)

```

